Multiple Numerical Assembly

Sections

Abstract

In this document, I want to provide a write-up of how this model is similar to and different from the versions previously implemented. This model is meant to provide assembly of multiple sites at the same time, that may, or may not, be connected in some fashion. Along the way, I will mention some of the input and output format that I am expect as documentation. At the end, I will be presenting some preliminary results. I especially want to use this as a vehicle to think about how to analyse those results.

Functions

library(dplyr)     # Data Manipulation
library(ggplot2)   # 2-D Plot
library(plotly)    # 3-D Plot
library(ggfortify) # used for biplots of PCAs
library(vegan)     # Ecological analysis mega-package
# library(RMTRCode2)

# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
    obj <- deparse(substitute(obj))
    if(exists(obj, envir = env)) {
        rm(list = obj, envir = env)
    }
}

Results, No Spatial Structure

First, we load up the preliminary results. In this case, we are loading a system in which 34 basal species and 66 consumer species form the pool for 10 unconnected environments. The pool and interaction matrix were assembled with the default parameters from Law and Morton’s 1996 work.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt-Result-Env10-None.RData")
)

Events

In total, 9320 events were used in these environments, with the species and environment invasion both randomly assigned. The number of arrival and extinction events were controlled to both be half of this number. We chose this number due to the coupon collecting problem. In particular, we use the result that the probability of encountering each species is bounded: \[\text{Pr}(\text{Draws} < n \log_{e} n + c n) \rightarrow \exp(-\exp(-c)) \text{ as } n \rightarrow \infty\] where \(n\) is the number of species and \(c\) is a constant. For our purposes, we choose \(c = 5\) so that we have a probability of about \(99.3\%\) of seeing each species in each environment. In practice, we failed to observe 2 species-environment combinations. Notably, nearly every species had at least one successful invasion; 72 did not.

The initial abundance was set to be 4000 times the elimination threshold, in line with work on minimum viable populations (Traill et al. 2007). The elimination threshold is admittedly more arbitrary, since it sets an effective individual-area relationship. For this calculation, I set it to \(10^{-4}\), in line with our previous calculations. This is large enough to avoid numerical difficulties from precision, while being low enough to represent a decent sized region.

Since each population is assembling simultaneously, I chose to use exponential waiting times for the inter-species arrival and extinction times. Note that these rates are shared between species and environments, but arrival and extinction are fully independent of each other. Species and environment affected in each event was chosen uniformly at random. The question then is how to set the rate.

To set the rate in this case, I chose to set it to the largest eigenvalue magnitude of the per-environment interaction matrices. This magnitude corresponds to the strongest response we can see from the interaction matrix and, if the interaction matrix is a good approximation for the Jacobian around a stable fixed point (which is not guaranteed), indicates the characteristic time scale of the decay to equilibrium. Hence, (overall) arrival rates and (overall) extinction rates should happen on the same timescale as the (largest) dynamics in the system. Since there are 10 environments, we should then expect that 10 characteristic time scales, on average, should occur in between arrival events in the same environment.

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[, c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Every vertical line is a species introduction or extinction by neutral dynamics.

Alpha Diversity plots

Intro

Perhaps more intriguing might be some sense of the biodiversity that we have in each system. We break the abundance results up by environment, then calculate the number of non-zero abundance curves at each time point. We also calculate the Shannon entropy (reminder: higher entropy means more uncertainty which means a flatter distribution).

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Richness = mean(Richness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

So richness hovers around a similar regime throughout the majority of the simulation. Note in this plot that we have emphasised one environmental curve and superimposed the mean in black. We do manage to reach heights of 11 species in one environment, but these heights are shortlived. Instead, we seem to observe a (time and environment averaged) value of 8.3677307. (If we consider the first 10,000 time units as burn-in, we instead see a value of 8.2536852.)

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Entropy has a similar, but highly erratic, behaviour. If one tries to follow one of the entropy curves, then one sees that they have fairly substantial periods of almost smooth behaviour followed by suddenly very noisy behaviour, and noise seems to be the dominant mode if one tries to examine the low alpha environment in the background. There are some easy to make predictions. Extinctions reduce the entropy in the system, as you become more certain about what remains. Analogously, arrivals increase the entropy. We can probably better see the relationship beyond these principles by plotting entropy against richness and connecting observations by environment and time.

Entropy-Richness

# ggplot2::ggplot(
#   Diversity %>% dplyr::filter(Environment < 4), 
#   ggplot2::aes(
#     x = Richness,
#     y = Entropy,
#     group = Environment,
#     color = Time
#   )
# ) + ggplot2::geom_path(
# ) + ggplot2::guides(
#   alpha = "none"
# )

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

It seems quite hard to tell, but there does not appear to be any particular orientation (clockwise, counter-clockwise) or similar pattern here.

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Evenness helps highlight that this is a high variance process but with a relatively constrained mean.

Environment Diversity

We can, of course, flip the idea on its head. Instead of examining the diversity of species within environments, we can look at the diversity of environments occupied by species. Since very few species end up occupying environments, we just look at richness. Unfortunately, this is quite a memory exhaustive task.

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

One immediately interesting trend here is that very few species are present across more than 5 environments at a given time. Indeed, only 2, 4, 5, 12, 13, 14, 22, 24, 29, 30, 32, 33, 35, 37, 43, 44, 47, 55, 57, 61, 68, 69, 71, 72, 76, 77, 78, 86 are ever present in more than 5 environments at once. We can also examine how long these periods occur for by species by tabulation. We block times so that entries are all of the same unit length, and round the average richness during the given time unit.

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

Beta Diversity

Principal Components Analysis

(For the desktop, the amount of data we generated is too much, so we need to reduce the amount of data we use. To do so we average over time blocks, here of length 100.)

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

We can perform a PCA and see if are data can be summarised by a small number of dimensions. As we shall see, constraining to the first 25 principal components does not harm the system much.

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

There is not actually a lot of dependence within the system it appears. Consider the amount of variance explained by the first six principal components (ordered, as is tradition, by amount of variation explained).

head(summary(PCA)$importance[2, ])

The amount of explained variation is as follows.

sum(summary(PCA)$importance[2, ])

The traditional biplot follows, but despite the seeming presence of patterns, so little of the variance is explained that is probably not worth further examination.

ggplot2::autoplot(PCA, loadings = TRUE, 
                  data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

We can try again with presence-absence data instead.

AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::mutate(
  dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCAPA)$importance[2, ])

So it is a bit better, but not by much.

ggplot2::autoplot(PCAPA, loadings = TRUE, 
                  data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

I should note that this is not a failure of the method persay. What this says to me is that dimension reduction of the system cannot reduce the system to a human-readable set of descriptors. The system is still being reduced (from 100 Species * 10 Environments to 25 Components to cover

sum(summary(PCA)$importance[2, ])

of the variation). My initial hypothesis for when the systems are coupled is that, as coupling increases, the number of principal components should decrease, since one system’s changes will be better predictors of the next system’s changes. (An interesting related question; do the number of principal components correlate with the amount of biodiversity in the system?

ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)

Beta Diversity Over Time

Since trying to reduce the system dimensionally does not work (well enough) a next attempt might be to consider how the pair-wise beta diversity changes over the course of the simulation. Initial problems include that there are quite a few measures of beta diversity as well as that the (square of the) number of environments will determine how many entries we need to consider.

To try this, we are going to reorganise our data into data frame with attributes of environment, time, and traits/species.

Environments <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    
    retval <- data.frame(Environment = toString(i),
                         Time = time
    )
    
    retval <- cbind(retval, env)
    colnames(retval) <- c("Environment", "Time", 
                          paste0("Basal", 1:34), 
                          paste0("Consumer", 35:100)
    )
    return(retval)
                          
  },
  abund = result$Abundance %>% data.frame(
  ) %>% dplyr::mutate(
    time = floor(time/100)*100
  ) %>% dplyr::group_by(
    time
  ) %>% dplyr::summarise(
    dplyr::across(.fns = ~ mean(.x))
  ),
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
  dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)

Cluster Analysis Over Time

One perhaps interesting idea is to try to group the environments by cluster by considering the abundances (or presence-absence) of the species present as traits. The question then is whether the environments have a tendency to attract to specific points or if they wander around each other without relation. (The latter is more neutral.) If they appear to be convergent (which so far would seem to disagree mostly with the alpha diversity analyses), then that would imply that dynamics determine the majority of the system’s fate, while if they instead wander more randomly, then they would appear to be dominated by the neutral mechanisms. (Of course, this is affected by the rate of the neutral mechanisms, so it exists on a definite continuum.) (Note, of course, that cluster analysis usually includes something like PCA, so we would not necessarily expect a different result.)

My working hypothesis for how to use the distance measures is that we need a metric measure (Kindt and Coe, 2005) with any identification scrubbed off (site, time collected). Kindt and Coe suggest that taking the square root of Bray-Curtis makes it metric and thus usable for mathematical distance analyses (such as hierarchical clustering). The vegan package recomends using Jaccard instead of Bray-Curtis, for the same reason.

EnvironmentDistance <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard")

EnvDistClust <- hclust(EnvironmentDistance)

plot(EnvDistClust, labels = FALSE)

Note that are major breakdowns extremely early yielding a large number of clusters but there are expected to be 10 environments (if historical contingency mattered most) or a small number of clusters (if they are converging). The large number of clusters says to me that the neutral dynamics are jumping rapidly from cluster to cluster.

We can try again with presence absence instead.

EnvironmentDistancePA <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)

EnvDistClustPA <- hclust(EnvironmentDistancePA)

plot(EnvDistClustPA, labels = FALSE)

Results, Line Spatial Structure

We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt-Result-Env10-Line.RData")
)

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[, c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Every vertical line is a species introduction or extinction by neutral dynamics.

Alpha Diversity plots

Intro

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Richness = mean(Richness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 1109 row(s) containing missing values (geom_path).
Warning: Removed 104 row(s) containing missing values (geom_path).

Entropy-Richness

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 1440 row(s) containing missing values (geom_path).
Warning: Removed 136 row(s) containing missing values (geom_path).

Environment Diversity

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

       Richness
Species     0     1     2     3     4     5     6     7     8     9    10
    1   80793     0     0     0     0     0     0     0     0     0     0
    2   80623     0     0     0     0     0     0     1     0     0   169
    3   80793     0     0     0     0     0     0     0     0     0     0
    4   80047     0     0     0     0     0     0     0     2     0   744
    5   79899     0     0     0     0     0     1     1     0     0   892
    6   80793     0     0     0     0     0     0     0     0     0     0
    7   80793     0     0     0     0     0     0     0     0     0     0
    8   80793     0     0     0     0     0     0     0     0     0     0
    9   80793     0     0     0     0     0     0     0     0     0     0
    10  80793     0     0     0     0     0     0     0     0     0     0
    11  80793     0     0     0     0     0     0     0     0     0     0
    12    680     0     0     0     0     0     0     0     0     0 80113
    13  80638     0     0     0     0     0     1     0     0     0   154
    14     84     0     0     0     0     0     0     0     0     0 80709
    15  80793     0     0     0     0     0     0     0     0     0     0
    16  80793     0     0     0     0     0     0     0     0     0     0
    17  80793     0     0     0     0     0     0     0     0     0     0
    18  80793     0     0     0     0     0     0     0     0     0     0
    19  80793     0     0     0     0     0     0     0     0     0     0
    20  80793     0     0     0     0     0     0     0     0     0     0
    21  80793     0     0     0     0     0     0     0     0     0     0
    22     62     0     0     0     0     0     0     0     0     0 80731
    23  80793     0     0     0     0     0     0     0     0     0     0
    24  80545     0     0     0     0     1     0     0     0     0   247
    25  80793     0     0     0     0     0     0     0     0     0     0
    26  80793     0     0     0     0     0     0     0     0     0     0
    27  80793     0     0     0     0     0     0     0     0     0     0
    28  80793     0     0     0     0     0     0     0     0     0     0
    29  80720     1     0     0     0     0     0     0     0     0    72
    30  69383     0     1     1     1     0     0     2     1     2 11402
    31  80793     0     0     0     0     0     0     0     0     0     0
    32    500     0     0     0     0     0     0     0     0     1 80292
    33  66478     0     0     1     3     4     1     1     0     2 14303
    34  80793     0     0     0     0     0     0     0     0     0     0
    35  80331     0     0     0     0     0     0     0     2     0   460
    36  80793     0     0     0     0     0     0     0     0     0     0
    37  46276     0     0     2     0     0     0     0     0     1 34514
    38  80793     0     0     0     0     0     0     0     0     0     0
    39  80793     0     0     0     0     0     0     0     0     0     0
    40  80793     0     0     0     0     0     0     0     0     0     0
    41  80793     0     0     0     0     0     0     0     0     0     0
    42  80793     0     0     0     0     0     0     0     0     0     0
    43  27626     0     0     1     0     2     0     1     1     0 53162
    44  60221     0     0     0     0     3     1     0     0     3 20565
    45  80793     0     0     0     0     0     0     0     0     0     0
    46  80793     0     0     0     0     0     0     0     0     0     0
    47  52181     0     0     1     0     1     1     0     0     0 28609
    48  80793     0     0     0     0     0     0     0     0     0     0
    49  80793     0     0     0     0     0     0     0     0     0     0
    50  80793     0     0     0     0     0     0     0     0     0     0
    51  80793     0     0     0     0     0     0     0     0     0     0
    52  80793     0     0     0     0     0     0     0     0     0     0
    53  80793     0     0     0     0     0     0     0     0     0     0
    54  80793     0     0     0     0     0     0     0     0     0     0
    55  78733     0     1     1     1     0     0     0     1     1  2055
    56  80793     0     0     0     0     0     0     0     0     0     0
    57  79019     0     0     0     0     0     0     0     0     0  1774
    58  80793     0     0     0     0     0     0     0     0     0     0
    59  80793     0     0     0     0     0     0     0     0     0     0
    60  80793     0     0     0     0     0     0     0     0     0     0
    61  80465     0     0     0     0     0     0     1     0     0   327
    62  80793     0     0     0     0     0     0     0     0     0     0
    63  80793     0     0     0     0     0     0     0     0     0     0
    64  80793     0     0     0     0     0     0     0     0     0     0
    65  80793     0     0     0     0     0     0     0     0     0     0
    66  80793     0     0     0     0     0     0     0     0     0     0
    67  80793     0     0     0     0     0     0     0     0     0     0
    68  79018     0     1     0     1     0     0     1     2     4  1766
    69  78129     0     1     0     0     2     0     2     1     0  2658
    70  80793     0     0     0     0     0     0     0     0     0     0
    71   1394     0     0     0     0     0     0     0     0     0 79399
    72    292     0     0     1     0     0     0     0     0     0 80500
    73  80793     0     0     0     0     0     0     0     0     0     0
    74  80793     0     0     0     0     0     0     0     0     0     0
    75  80793     0     0     0     0     0     0     0     0     0     0
    76  78672     1     0     0     0     1     0     1     0     1  2117
    77  79764     0     0     0     0     0     0     0     0     0  1029
    78  79252     0     1     0     1     0     0     0     0     0  1539
    79  80793     0     0     0     0     0     0     0     0     0     0
    80  80793     0     0     0     0     0     0     0     0     0     0
    81  80793     0     0     0     0     0     0     0     0     0     0
    82  80793     0     0     0     0     0     0     0     0     0     0
    83  80793     0     0     0     0     0     0     0     0     0     0
    84  80793     0     0     0     0     0     0     0     0     0     0
    85  80793     0     0     0     0     0     0     0     0     0     0
    86  64857     0     0     1     1     3     0     1     1     1 15928
    87  80793     0     0     0     0     0     0     0     0     0     0
    88  80793     0     0     0     0     0     0     0     0     0     0
    89  80793     0     0     0     0     0     0     0     0     0     0
    90  80793     0     0     0     0     0     0     0     0     0     0
 [ reached getOption("max.print") -- omitted 10 rows ]

Beta Diversity

Principal Components Analysis

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.19145 0.11089 0.10844 0.07943 0.05993 0.04530 

sum(summary(PCA)$importance[2, ])

[1] 0.99996

ggplot2::autoplot(PCA, loadings = TRUE, 
                  data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::mutate(
  dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCAPA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.21745 0.12196 0.08092 0.05942 0.05217 0.04471 

ggplot2::autoplot(PCAPA, loadings = TRUE, 
                  data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

sum(summary(PCA)$importance[2, ])

[1] 0.99996

ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)

Beta Diversity Over Time

Environments <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    
    retval <- data.frame(Environment = toString(i),
                         Time = time
    )
    
    retval <- cbind(retval, env)
    colnames(retval) <- c("Environment", "Time", 
                          paste0("Basal", 1:34), 
                          paste0("Consumer", 35:100)
    )
    return(retval)
                          
  },
  abund = result$Abundance %>% data.frame(
  ) %>% dplyr::mutate(
    time = floor(time/100)*100
  ) %>% dplyr::group_by(
    time
  ) %>% dplyr::summarise(
    dplyr::across(.fns = ~ mean(.x))
  ),
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
  dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)

Cluster Analysis Over Time

EnvironmentDistance <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard")

EnvDistClust <- hclust(EnvironmentDistance)

plot(EnvDistClust, labels = FALSE)

EnvironmentDistancePA <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)

EnvDistClustPA <- hclust(EnvironmentDistancePA)

plot(EnvDistClustPA, labels = FALSE)